Conduction-band tight-binding description for Si applied to P donors 
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A tight-binding parametrization for silicon, optimized to correctly reproduce effective masses as 
well as the reciprocal space positions of the conduction-band minima, is presented. The reliability 
of the proposed parametrization is assessed by performing systematic comparisons between the 
descriptions of donor impurities in Si using this parametrization and previously reported ones. The 
spectral decomposition of the donor wavefunction demonstrates the importance of incorporating full 
band effects for a reliable representation, and that an incomplete real space description results from 
a truncated reciprocal space expansion as proposed within the effective mass theory. 



Advances in semiconductor device fabrication, partic- 
ularly Si-based devices, have benefited from the progres- 
sive miniaturization and integration of their constituent 
parts, including detailed control of the doping process. 
The theoretical approach utilized for device modeling 
should be able to resolve charge-density variations on 
an atomic scale, and in this respect tight-binding (TB) 
models are particularly adequate, as they provide atom- 
istic descriptions of structural and electronic properties of 
solids. 1 Since the pioneering work of Slater and Koster,— 
TB is conceived as an empirical method, where orbital 
energies and hoppings are parameters to be adjusted in 
order to reproduce relevant properties in the band struc- 
ture of the solid. TB methods rely on the choice of a suit- 
able parametrization of the Hamiltonian, which includes 
the choice of a basis set, range of hopping coupling, etc. 
For the group IV and III-V semiconductors one might 
expect that the minimal sp 3 basis with first-neighbors 
coupling should suffice to describe the essential physics. 
However, this model is not able to reproduce the gaps 
of the indirect gap materials like Si and AlAs.— A simple 
scheme to correct this deficiency consists in introducing 
an excited s-like orbital, s* , to improve the description of 
the conduction band (CB). 4 Another limitation of first- 
neighbor models based on s and p orbitals alone is that 
they predict an infinite value of the transverse effective 
mass at the X point. This anomaly may be overcome 
either by adding the five d orbitals in the first-neighbor 
basis set,— £ or by inclusion of second- neighbor (2nn) in- 
teractions in the sp 3 s* model.— 

Comparison between the different TB parametriza- 
tions for Si available in the literature shows that the one 
that best reproduces the relevant bulk material proper- 
ties is the one proposed by Klimeck et al^ based on the 
first and second-neighbors sp 3 s* model. The reciprocal 
space position of the CB minima, however, is not well 
fitted. Silicon is an indirect gap material, with 6 CB 
minima along the equivalent A lines, 85% of the way be- 
tween r and X . In Ref. 0, a genetic algorithm (GA) 
fitting was adopted, and the target position for the min- 



ima was considered as 75% between T and X. For some 
applications^ it is known that a shift in the position of 
the CB minimum may lead to unsatisfactory results. We 
present here an improved parametrization, also obtained 
through the GA methodology^ reproducing the correct 
positions of the CB edges. 

GA-based optimization employs stochastic methods 
which do not require constrains on continuity of solu- 
tion space. In the case considered here, the TB parame- 
ters are varied in order to get an optimal set that repro- 
duces given material properties, denoted as target values. 
Details on the method are discussed in Ref. We re- 
fer to the parametrization published there as P075, and 
to the one proposed here as P085. Both parametriza- 
tions give the fc-space positions of the six band minima 
at six equivalent points along A lines, with the position 
of the minimum at A m i„ = 0.75(27r/asi) for P075 and 
A mm = 0.85(27r/a si ) for P085, where a si = 5.43lA is 
the conventional cubic lattice parameter for Si. Table U 
presents the TB parameters for Si for both parametriza- 
tions. The main difference between them consists in al- 
lowing non-zero 2nn hoppings V S s*(110) and 14* s *(110) 
in the optimization set for P085. All 2nn hoppings 
are consistently smaller than the first-neighbors. The 
2nn hopping parameters were determined to adjust the 
finer details of the target properties, and were not con- 
strained to have values or signs expected from physical 
considerations^ 

The input and the calculated properties are presented 
in Table |TT1 where the first column shows the mate- 
rial properties to be represented by the TB model and 
the second column shows the corresponding experimen- 
tal values^ These constitute the input targets in the GA 
code. The remaining columns give the calculated proper- 
ties and the respective deviations from the experimental 
targets for P075 and P085, as well as for a recently pro- 
posed first-neighbor sp 3 d 5 s* parametrization, which we 
denote by Plnn£ We note from Table |TT] that the P075 
and P085 parametrizations give consistently better agree- 
ment with the target values than the Plnn parametriza- 
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TABLE I: Parameters for the TB models (in eV). 



Parameter 


P075 


P085 


£ s (000) 


-4.81341 


-4.848054 


£p(000) 


1.77563 


1.787118 


E s . (000) 


6.61342 


5.608014 




-8.33255 


-8.259704 


V xx (ill) 


1.69916 


1.697556 




5.29091 


5.351079 




5.86140 


5.822197 


V s *p (222) 


4.88308 


4.864480 


Aso 


0.04503 


0.014905 


V„(110) 


0.01591 


0.029958 


VW(110) 


0.00000 


0.191517 


T4 S .(110) 


0.00000 


0.007036 


V ax {110) 


0.08002 


0.161749 


14,(011) 


1.31699 


0.885988 


K*a:(110) 


-0.00579 


-0.095653 


14-4011) 


0.50103 


0.966257 


14.(110) 


0.00762 


0.037296 


14,(011) 


-0.10662 


-0.132810 


T4„(110) 


0.55067 


0.619876 


K„(011) 


-2.27784 


-2.496288 



TABLE II: Optimization targets and optimized material 
properties for the P075, P085 and Plnn models. Except for 
A m i„, which is specific for different models, the target values 
correspond to experimental data given in Ref. Hbl . 



Property 


Target 


P075 


%dev 


P085 


%dev 


Plnn %dev 




0.750 


0.758 


1.067 












0.850 






0.8480 


-0.235 


0.813 


-4.35 


El 


3.350 


3.353 


0.089 


3.350 


-0.013 


3.399 


1.44 




1.130 


1.129 


-0.050 


1.130 


-0.042 


1.131 


0.09 


m*xi 


0.916 


0.916 


-0.030 


0.916 


0.050 


0.891 


-2.73 




0.191 


0.191 


0.020 


0.191 


0.007 


0.201 


-5.23 


rof h [001] 


-0.204 


-0.198 


3.082 


-0.204 


-0.060 


-0.214 


-4.90 


KU011] 


-0.147 


-0.146 


0.525 


-0.148 


-0.568 


-0.152 


-3.40 


n»H,[Hl] 


-0.139 


-0.139 


0.395 


-0.140 


-0.610 


-0.144 


-3.60 


m^[001] 


-0.275 


-0.285 


-3.643 


-0.277 


-0.786 


-0.276 


-0.36 


m* hh [011] 


-0.579 


-0.581 


-0.338 


-0.574 


0.869 


-0.581 


-0.34 


m* hh [in] 


-0.738 


-0.737 


0.119 


-0.727 


1.466 


-0.734 


0.54 


m* so 


-0.234 


-0.237 


-1.487 


-0.239 


-2.162 


-0.246 


-5.13 


A so 


0.015 


0.145 


-0.067 


0.015 


0.030 


0.016 


4.90 



tion. This statement cannot be generalized as a fun- 
damental advantage of the second-neighbor sp 3 s* model 
over the nearest-neighbor sp 3 d 5 s* model. The Plnn sei£ 
satisfies additional requirements beyond bulk behavior 
properties. On-going work indicates that the sp 3 d 5 s* 
model can be fit to match the Si bulk properties just 
as well as the second-neighbor sp 3 s* model. The pri- 
mary advantage of the nearest-neighbor sp 3 d 5 s* model is 



its straightforward incorporation of strain distortions. 11 
In the donor description below, no strain distortions are 
considered, and the second-neighbor model P085 pro- 
vides an unprecedented representation of the Si CB. 

The reliability of TB sp 3 s* second-neighbors 
parametrizations has been recently verified by studies of 
shallow donors in GaAsi^ and in SiJ£ We perform the 
same kind of study here, and discuss the effect of the 
CB minimum position in different aspects of the donor 
problem. We write the Hamiltonian for the impurity 
problem as 14 H = £ J2 pu h^c]^ + J2 t .„ U{Ri)c\ v c iv 

where i and j label the atomic sites, /1 and v denote 
the atomic orbitals and Ri is the distance between site 
i and the impurity site. The impurity potential is taken 



as a screened Coulomb potential, U(Ri) 



(e = 12.1 for Si). At the impurity site it is assigned the 
value U(Ri = 0) = — Uq, a parameter describing central 
cell effects characteristic of the substitutional species, 
and taken here as an adjustable parameter. We do not 
include spin-orbit corrections in the present calculations. 

The eigenstates of H are determined for a system 
where a single impurity is placed in a cubic supercell 
containing N — 8L 3 atoms in the diamond structure, 
where L is the length of the supercell edge in units of 
asi. We adopt periodic boundary conditions, and large 
supercells 1 ^ (up to 10 6 atoms) were treated within a 
variational schem o 15 ' 16 where the ground state wavefunc- 
tion and binding energy El are obtained by minimizing 

(H — s re f) 2 The reference energy £ re / is chosen 

within the gap, nearest to the CB edge. 

The eigenfunctions of H in the basis of atomic-like or- 
bitals are written as |*Ts(r)> = Z) w a ^l^( r ~ 
where the expansion coefficients di„ give the probabil- 
ity amplitude of finding the electron in the orbital v at 
site The overall charge distribution is conveniently 
described through the TB envelope function squared?^ 



(i) 



In Fig. [TJ a convergence study of the donor ground 
state binding energy as a function of the supercell size 
L is presented for P075 and P085. One can observe 
in Fig. [IJa) that for supercell sizes L > 25 the calcu- 
lated binding energies reproduce the experimental value 
(E exp = 45.6 meV) taking U = 1.48 and U = 1.26 eV 
for P075 and P085, respectively. We denote these values 
by Up, as they are determined by tuning the on-site po- 
tential — Uq in order to give the converged value of E u in 
agreement with experiment for P donors in Si. Fig. |Tfb) 
presents a plot of ln[(i?£ — E exp ) / E exp ] vs L. The linear 
behavior obtained here indicates the same exponential 
convergence of El to the experimental value for both 
parametrizations: El ~ E exp + Ee~ L / x . 

We determine the first excited state energy^ by vary- 



ing the value of e re f in ( 



(H - SrefY 



^ ) . Keeping 



Uq = 1-26 eV, we obtain the binding energy of the first 
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FIG. 1: (Color online) (a) Convergence of the donor ground 
state binding energy towards the experimental value (dotted 
line gives E exp = 45.6 meV) with the supercell size L, for 
U = Up, namely 1.48 eV and 1.26 eV for P075 and P085 
respectively, (b) Same data plotted as ln[(i5_L — E exp ) / E exp ] 
vs L. The linear behavior (the line is a best fit for all data 
points) indicates that the convergence of the binding energy 
with the supercell size L follows the same exponential law for 
both parametrizations. 



excited state to be 30.2 meV, in good agreement with the 
P075 result, 32.4 meV. 

The behavior of the binding energy with Uq is pre- 
sented in Fig. (Ha), where the dotted lines indicate the 
value of Up for each parametrization leading to Et = 
E exp . As noted in Ref. 1 131 in the weak perturbation limit 
the binding energies converge to the effective mass theory 
(EMT) prediction in its simplest formulation 1 ^ (single- 
valley approximation), ~ 30 meV. It is interesting that 
this behavior is the same for both parametrizations. As 
Uo increases, P085 tends to give higher binding energies 
than P075, resulting in a smaller value of Up for P085. 

One can also characterize the donor ground state by 
its orbital averaged spectral weight^ at k = A m j n , 

w(A mm ) = ^EiLiE^-^Woj 



plotted in 



Fig. [2jb) as a function of Uq. The EMT approach pre- 
sumes that W(A m m) ~ 1, allowing the donor state to 
be well described in a basis of Bloch states at the CB- 
edge k points. However, one can notice in Fig. [3] (b) 
that this is not the case: Even for the smallest values 
of Uq the spectral weights at A mirl are well below satu- 
ration (one) for both parametrizations, implying that an 
incomplete description may result from EMT in this case. 
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FIG. 2: (Color online) (a) Binding energy of the ground state 
as a function of the on-site perturbation strength Uq. The 
dotted lines indicate the value Uo = Up that reproduces the 
experimental Si:P A\ state binding energy, (b) Total spectral 
weight at the CB edges for the ground impurity state. 
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FIG. 3: (Color online) TB envelope function squared for the 
donor ground state along the [110] direction. The lines are 
the corresponding effective mass \ip\ 2 results. 
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The spectral weight at U P is 0.32 for P075 and 0.30 for 
P085. These relatively low spectral weights indicate that 
multiple k points, other than those corresponding to the 
six CB-edges of Si, contribute to the donor wavefunction 
expansion within any reciprocal-space based approach. 

Within single-valley EMT the ground state for donors 
in Si is six-fold degenerated Valley-orbit interactions^ 
lead to a non-degenerate ground state wavefunction of 
A\ symmetry, 

1 6 

^(r) = - 7 =^F /1 (rK(r)e ik -- r , (2) 

where </>u(r) = w Al (r)e 4k,J ' r are the pertinent Bloch 
wavefunctions, and the envelope functions are given by 
F z (r) = (1/^) e -[(* 2 + y 2 )/a 2 +,> 2 ] 1/2 for ^ = z and 
equivalently for the other /x values. The effective Bohr 
radii for Si are a — 2.51 nm and b = 1.44 nm^ Fig. [3] 
presents a comparison between the TB envelope func- 
tion calculated from ([I]) along the [110] direction, for 
both P075 and P085 parametrizations (data points) with 
the corresponding EMT results obtained from @ (solid 
lines), with a and b are given above, but with a different 
normalization to conciliate the TB and EMT wavefunc- 
tions on the same scale. The values of used in ([TJ are 
consistent with the respective reciprocal space location of 
the CB minima. Note that the oscillatory behavior due 
to interference among the plane-wave parts of the six (j)^ 
is well captured by the TB envelope function. The pe- 
riod of the oscillations is different for P075 and P085, as 
given by the corresponding wavevectors. 

Good agreement between TB and EMT is restricted 
to distances from the impurity site larger than ~ 1 nm. 
This means that at large distances the EMT expansion 
of the donor wavefunction in only six k points, as given 



in Eq. (J2]) , is capable of reproducing its main features (ex- 
cept for normalization, of course). Closer to the impurity, 
particularly at the impurity site, the TB results become 
much larger than the EMT prediction. 

The wavefunctions obtained from P085 and P075 agree 
reasonably well in the central cell region. One way 
to quantify this agreement is through the probability 
to find the donor electron inside of a sphere of radius 
Re- Q(R C ) = Er<b. c MR)?- Taking for R c the 
2nn distance, the ratio of Q(R C ) obtained from the two 
parametrizations is Qpo85(-Rc)/Qpo75(-Rc) = 1-15. 

In summary, we find good agreement between the re- 
sults obtained within P085 and P075 for (i) the expo- 
nential convergence law for the ground state binding en- 
ergy with supercell size, (ii) the binding energy of the 
first excited state, (iii) the spectral weight of the ground 
state wavefunction at A m i„, (iv) the probability that the 
donor electron is within the central cell up to the im- 
purity's 2nn. Both parametrizations also capture the 
donor wavefunction oscillations predicted within EMT; 
the main difference regards the period of the oscillations. 
In applications where the quantitative aspects of the os- 
cillatory behavior of the wavefunction is important, the 
P085 parametrization is thus capable of providing a bet- 
ter description. The importance to represent the P impu- 
rity wavefunction in a full band, atomistic representation 
is demonstrated through its spectral decomposition, and 
by direct comparison between the TB results and EMT. 
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